//-----------------------Bibliotecas-----------------------------------------
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "Matrizes.h"
//---------------------------------------------------------------------------
#define MAX(x,y) ((x)>(y)?(x):(y))
#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//----------------Algoritmos de Manipulacao de Matrizes----------------------
//---------------------------------------------------------------------------
//--------------------------Corpo das Funcoes--------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
float *CriaVetor(int NumeroElementos)
{
   float *Vetor;

   Vetor = new float[NumeroElementos];
   if(Vetor==NULL)
      return(NULL);

   return(Vetor);
}
//---------------------------------------------------------------------------
float *ApagaVetor(float *A)
{
  if (A == NULL) return (NULL);
  delete[] A;
  return (NULL);
}
//----------------------------------------------------------------------------
float **CriaMatriz(int NumeroLinhas, int NumeroColunas)
{
   float **Matriz;

   Matriz = new float*[NumeroLinhas];

   if(Matriz==NULL) return(NULL);

   for(int i=0;i<NumeroLinhas;i++)
   {
      Matriz[i] = new float[NumeroColunas];
   }

   return(Matriz);
}
//---------------------------------------------------------------------------
float **ApagaMatriz(float **A,int NumeroLinhas, int NumeroColunas)
{
  if (A == NULL) return (NULL);
  for (int i=0; i<NumeroLinhas; i++)
     delete[] A[i];
  delete[] A;

  return (NULL);
}
//----------------------------------------------------------------------------
float **MultiplicaMatriz(float **A,int LinhasA,int ColunasA,float **B,int LinhasB,int ColunasB)
{
   float **Resultado;
   long double Apoio;

   Resultado = CriaMatriz(LinhasA,ColunasB);
   //ColunasA == LinhasB
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ColunasB;j++)
      {
         Apoio = 0.0;
         for(int k=0;k<LinhasB;k++)
            Apoio+=(long double)(A[i][k]*B[k][j]);
         Resultado[i][j] = (float) Apoio;
      }
   return(Resultado);
}
//---------------------------------------------------------------------------
float **MultiplicaConstante(float **A,int LinhasA,int ColunasA,float B)
{
   float **Resultado;

   Resultado = CriaMatriz(LinhasA,ColunasA);
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ColunasA;j++)
         Resultado[i][j]=B*A[i][j];
   return(Resultado);
}
//---------------------------------------------------------------------------
float  *MultiplicaConstante(float  *A,int ElementosA,float B)
{
   float *Resultado;

   Resultado = CriaVetor(ElementosA);
   for(int i=0;i<ElementosA;i++)
         Resultado[i] = B*A[i];
   return(Resultado);
}
//---------------------------------------------------------------------------
float **SomaMatriz(float **A,int LinhasA,int ColunasA,float **B)
{
   float **Resultado;

   Resultado = CriaMatriz(LinhasA,ColunasA);
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ColunasA;j++)
         Resultado[i][j]=A[i][j]+B[i][j];
   return(Resultado);
}
//---------------------------------------------------------------------------
float **SubtraiMatriz(float **A,int LinhasA,int ColunasA,float **B)
{
   float **Resultado;

   Resultado = CriaMatriz(LinhasA,ColunasA);
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ColunasA;j++)
         Resultado[i][j]=A[i][j]-B[i][j];
   return(Resultado);
}
//----------------------------------------------------------------------------
float CalculaSomaElementos(float **Matriz,int NumLinhas,int NumColunas)
{
   int i,j;
   float Resultado = 0;

   for(i=0;i<NumLinhas;i++)
      for(j=0;j<NumColunas;j++)
      {
         Resultado += Matriz[i][j];
      }
   return Resultado;
}
//---------------------------------------------------------------------------
float   CalculaSomaElementos(float *Vetor,int NumElementos)
{
   int i;
   float Resultado = 0;

   for(i=0;i<NumElementos;i++)
      Resultado += Vetor[i];

   return Resultado;
}
//---------------------------------------------------------------------------
float CalculaMediaElementos(float **Matriz,int NumLinhas,int NumColunas)
{
   float Retorno;
   int i,j;
   float Apoio;

   Retorno = 0.0f;
   for(i=0;i<NumLinhas;i++)
   {
      Apoio = 0.0f;
      for(j=0;j<NumColunas;j++)
         Apoio += Matriz[i][j];
      Retorno += Apoio/NumColunas;
   }
   return (Retorno/NumLinhas);
}
//---------------------------------------------------------------------------
float CalculaMediaElementos(float *Vetor,int NumElementos)
{
   float Retorno;
   int i;

   Retorno = 0.0f;
   for(i=0;i<NumElementos;i++)
      Retorno += Vetor[i];

   return Retorno;
}
//---------------------------------------------------------------------------
float CalculaMaximoElemento(float **Matriz,int NumLinhas,int NumColunas)
{
   int i,j;
   float Retorno;

   Retorno = Matriz[0][0];
   for(i=0;i<NumLinhas;i++)
      for(j=0;j<NumColunas;j++)
         if(Retorno<Matriz[i][j])
            Retorno = Matriz[i][j];
   return Retorno;
}
float   CalculaMaximoElemento(float *Vetor,int NumElementos)
{
   int i;
   float Retorno;

   Retorno = Vetor[0];
   for(i=0;i<NumElementos;i++)
         if(Retorno<Vetor[i])
            Retorno = Vetor[i];
   return Retorno;
}
//---------------------------------------------------------------------------
float   CalculaMinimoElemento(float **Matriz,int NumLinhas,int NumColunas)
{
   int i,j;
   float Retorno;

   Retorno = Matriz[0][0];
   for(i=0;i<NumLinhas;i++)
      for(j=0;j<NumColunas;j++)
         if(Retorno>Matriz[i][j])
            Retorno = Matriz[i][j];
   return Retorno;
}
float CalculaMinimoElemento(float *Vetor,int NumElementos)
{
   int i;
   float Retorno;

   Retorno = Vetor[0];
   for(i=0;i<NumElementos;i++)
         if(Retorno>Vetor[i])
            Retorno = Vetor[i];
   return Retorno;
}
//---------------------------------------------------------------------------
float   CalculaVarianciaElementos(float **Matriz,int NumLinhas,int NumColunas)
{
 float Media;
 float Variancia;
 int i,j;

 Media = CalculaMediaElementos(Matriz,NumLinhas,NumColunas);
 for(i=0;i<NumLinhas;i++)
     for(j=0;j<NumColunas;j++)
        Variancia += (Matriz[i][j]-Media)*(Matriz[i][j]-Media);
 Variancia = Variancia/(NumLinhas*NumColunas-1);
 return Variancia;
}
//---------------------------------------------------------------------------
float **CalculaTransposta(float **A,int LinhasA,int ColunasA)
{
   float **Resultado;

   Resultado = CriaMatriz(ColunasA,LinhasA);
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ColunasA;j++)
         Resultado[j][i]=A[i][j];
   return(Resultado);
}
//---------------------------------------------------------------------------
int FatoracaoLU(float **A,int LinhasMatriz,float **Permutacao)
{
   int i,j,k,l;
   float temp;
   int temp2;
   int Sinal = 1;

   if(Permutacao!=NULL)
   {
      for(i=0;i<LinhasMatriz;i++)
         for(j=0;j<LinhasMatriz;j++)
            if(i==j)
               Permutacao[i][j] = 1.0;
            else
               Permutacao[i][j] = 0.0;
   }


   for(i=0;i<LinhasMatriz-1;i++)//colunas
   {
   //-------------------------Pivoteamento Parcial---------------------------
       temp=A[i][i];
       temp2=0;
       for(l=i+1;l<LinhasMatriz;l++)
       {
          if(temp*temp<A[l][i]*A[l][i])
          {
             temp = A[l][i];
             temp2 = l;
          }
       }
       if(temp2!=0)
       {
          Sinal = -Sinal;
          if(Permutacao!=NULL)
          {
             for(l=0;l<LinhasMatriz;l++)
             {
               temp=Permutacao[i][l];
               Permutacao[i][l]=Permutacao[temp2][l];
               Permutacao[temp2][l]=temp;
             }
          }
          for(l=i;l<LinhasMatriz;l++)
          {
            temp=A[i][l];
            A[i][l]=A[temp2][l];
            A[temp2][l]=temp;
          }
       }
   //------------------------------------------------------------------------
      for(j=i+1;j<LinhasMatriz;j++)    //linhas da coluna i
      {
    	 if(A[j][i]==0&&A[i][i]==0){
    	 temp =0;
    	 }else{
    	 temp = (A[j][i]/A[i][i]);
    	 }
         for(k=i+1;k<LinhasMatriz;k++)  //colunas da linha j
         {
            A[j][k]-=(A[i][k]*temp);
         }
         A[j][i]=temp;
      }
   }
   return(Sinal);
}
//---------------------------------------------------------------------------
float **RemoveColuna(float **Matriz,int NumLinhas, int NumColunas,int ColunaRemovida)
{
   float **Resultado;
   float temp;
   int i,j,k;

   Resultado = CriaMatriz(NumLinhas,NumColunas-1);
   k = 0;
   for(i=0;i<NumColunas;i++)
      if(i!=ColunaRemovida)
      {
         for(j=0;j<NumLinhas;j++)
         {
            Resultado[j][k] = Matriz[j][i];
         }
         k++;
      }
   return(Resultado);
}
//---------------------------------------------------------------------------
float **RemoveLinha(float **Matriz,int NumLinhas, int NumColunas,int LinhaRemovida)
{
   float **Resultado;
   float temp;
   int i,j,k;

   Resultado = CriaMatriz(NumLinhas-1,NumColunas);
   k = 0;
   for(i=0;i<NumLinhas;i++)
      if(i!=LinhaRemovida)
      {
         for(j=0;j<NumColunas;j++)
         {
            Resultado[k][j] = Matriz[i][j];
         }
         k++;
      }
   return(Resultado);
}
//---------------------------------------------------------------------------
float **CalculaCofatores(float **Matriz,int NumLinhas)
{
   float **Cofatores,
         **MatrizTemp,
         **MatrizTemp2;
   float Det;
   int i,j;

   Cofatores = CriaMatriz(NumLinhas,NumLinhas);

   for(i=0;i<NumLinhas;i++)
   {
	   MatrizTemp  = RemoveLinha(Matriz,NumLinhas,NumLinhas,i);
	   for(j=0;j<NumLinhas;j++)
	   {
		   MatrizTemp2 = RemoveColuna(MatrizTemp,NumLinhas-1,NumLinhas,j);
		   Det = CalculaDeterminante(MatrizTemp2,NumLinhas-1);
		   MatrizTemp2 = ApagaMatriz(MatrizTemp2,NumLinhas-1,NumLinhas-1);
		   if((i+j)%2){
			   Cofatores[i][j] = -Det;
		   }
		   else{
			   Cofatores[i][j] =  Det;
		   }
	   }
	   MatrizTemp  = ApagaMatriz(MatrizTemp,NumLinhas-1,NumLinhas);

   }
   return(Cofatores);
}
//---------------------------------------------------------------------------
float **CalculaAdjunta(float **A,int NumLinhas)
{
   float **Cofatores;
   float **Adjunta;

   Cofatores = CalculaCofatores(A,NumLinhas);
   Adjunta   = CalculaTransposta(Cofatores,NumLinhas,NumLinhas);

   ApagaMatriz(Cofatores,NumLinhas,NumLinhas);

   return(Adjunta);
}
//---------------------------------------------------------------------------
float *SolucaoSistemaLU(float **MatrizFatoradaLU,int NumLinhas,float **Permutacao,float *Vetor)
{
   return(NULL);
}

//---------------------------------------------------------------------------

float **CalculaInversa(float **Matriz,int NumLinhas)
{
   float **Inversa,
         **Adjunta;
   float InvDet;

   Adjunta = CalculaAdjunta(Matriz,NumLinhas);
   InvDet  = 1/(CalculaDeterminante(Matriz,NumLinhas));
   Inversa = MultiplicaConstante(Adjunta,NumLinhas,NumLinhas,InvDet);

   Adjunta = ApagaMatriz(Adjunta,NumLinhas,NumLinhas);

   return(Inversa);
}
//---------------------------------------------------------------------------
void EliminacaoGaussJordan(float **A,float **B,int LinhasMatriz,int ColunasMatriz)
{
   int i,j,k,l;
   float temp;
   int temp2;

   //Aplica eliminacao de Gauss em <A>*<X>=<B> , onde X=Vetor[0][LinhasMatriz]
   for(i=0;i<ColunasMatriz-1;i++)//colunas
   {

   //-------------------------Pivoteamento Parcial---------------------------
       temp=A[i][i];
       temp2=0;
       for(l=i+1;l<LinhasMatriz;l++)
       {
          if(temp*temp<A[l][i]*A[l][i])
          {
             temp  = A[l][i];
             temp2 = l;
          }
       }
       if(temp2!=0)
       {
          for(l=i;l<ColunasMatriz;l++)
          {
            temp=A[i][l];
            A[i][l]=A[temp2][l];
            A[temp2][l]=temp;
          }
          if(B!=NULL)
          {
             temp=B[i][0];
             B[i][0]=B[temp2][0];
             B[temp2][0]=temp;
          }
       }
   //------------------------------------------------------------------------

      for(j=i+1;j<LinhasMatriz;j++)    //linhas da coluna i
      {
         for(k=i+1;k<ColunasMatriz;k++)  //colunas da linha j
         {
            A[j][k]-=(A[i][k]*(A[j][i]/A[i][i]));
         }
         if(B!=NULL)
            B[j][0]-=(B[i][0]*(A[j][i]/A[i][i]));
         A[j][i]=0.0f;
      }

   }
}
//---------------------------------------------------------------------------

//ISSO É UM GATO!
float **matrizMenor(float **m, float **m2, int c, int n){
        int i, j;
        float **p=m2, *p2;
        for(i = 1; i < n; i++){
                p2 = *p++;
                for(j = 0; j < n; j++)
                        if(j==c) continue;
                        else *p2++ = m[i][j];
        }
        return m2;
}
//Retorna o determinante da matriz m de ordem n por Laplace
int det(float **m, int n){
        //Caso fundamental
        if(n == 1)
                return **m;

        //Criando a matriz do cofator
        int c,d=0;
        float cof;
        float **m2 = CriaMatriz(n,n);

        //Para cada coluna
        for(c=0; c<n; c++){
                //Calculando o cofator
                cof = (c%2?-1:1) * det(matrizMenor(m, m2, c, n), n-1);
                //Somando a parcela ao determinante final
                d += m[0][c]*cof;
        }

        //Liberando o espaço alocado da matriz do cofator
        ApagaMatriz(m2, n-1,n-1);

        //Retornando o determinante
        return d;
}



float CalculaDeterminante(float **Matriz,int NumLinhas)
{
	/*
   int Sinal,
       i,j;

   float **MatrizFatoradaLU;
   float Determinante = 1;

   MatrizFatoradaLU = CriaMatriz(NumLinhas,NumLinhas);
   for(i=0;i<NumLinhas;i++)
      for(j=0;j<NumLinhas;j++)
      {
         MatrizFatoradaLU[i][j] = Matriz[i][j];
      }

   Sinal = FatoracaoLU(MatrizFatoradaLU,NumLinhas,NULL);

   for(i=0;i<NumLinhas;i++){
	   Determinante = (Determinante)*(MatrizFatoradaLU[i][i]);
   }
   Determinante = Determinante*Sinal;

   MatrizFatoradaLU = ApagaMatriz(MatrizFatoradaLU,NumLinhas,NumLinhas);
*/


   return det(Matriz,NumLinhas);
}
//---------------------------------------------------------------------------
/*
 * svdcomp - SVD decomposition routine.
 * Takes an mxn matrix a and decomposes it into udv, where u,v are
 * left and right orthogonal transformation matrices, and d is a
 * diagonal matrix of singular values.
 *
 * This routine is adapted from svdecomp.c in XLISP-STAT 2.1 which is
 * code from Numerical Recipes adapted by Luke Tierney and David Betz.
 *
 * Input to dsvd is as follows:
 *   a = mxn matrix to be decomposed, gets overwritten with u
 *   m = row dimension of a
 *   n = column dimension of a
 *   w = returns the vector of singular values of a
 *   v = returns the right orthogonal transformation matrix
*/

static double PYTHAG(double a, double b)
{
    double at = fabs(a), bt = fabs(b), ct, result;

    if (at > bt)       { ct = bt / at; result = at * sqrt(1.0 + ct * ct); }
    else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); }
    else result = 0.0;
    return(result);
}


int DecomposicaoSVD(float **a, int m, int n, float *w, float **v)
{
    int flag, i, its, j, jj, k, l, nm;
    double c, f, h, s, x, y, z;
    double anorm = 0.0, g = 0.0, scale = 0.0;
    double *rv1;

    if (m < n)
    {
        fprintf(stderr, "#rows must be > #cols \n");
        return(0);
    }
    rv1 = (double *)malloc((unsigned int) n*sizeof(double));

    for (i = 0; i < n; i++)
    {
        l = i + 1;
        rv1[i] = scale * g;
        g = s = scale = 0.0;
        if (i < m)
        {
            for (k = i; k < m; k++)
                scale += fabs((double)a[k][i]);
            if (scale)
            {
                for (k = i; k < m; k++)
                {
                    a[k][i] = (float)((double)a[k][i]/scale);
                    s += ((double)a[k][i] * (double)a[k][i]);
                }
                f = (double)a[i][i];
                g = -SIGN(sqrt(s), f);
                h = f * g - s;
                a[i][i] = (float)(f - g);
                if (i != n - 1)
                {
                    for (j = l; j < n; j++)
                    {
                        for (s = 0.0, k = i; k < m; k++)
                            s += ((double)a[k][i] * (double)a[k][j]);
                        f = s / h;
                        for (k = i; k < m; k++)
                            a[k][j] += (float)(f * (double)a[k][i]);
                    }
                }
                for (k = i; k < m; k++)
                    a[k][i] = (float)((double)a[k][i]*scale);
            }
        }
        w[i] = (float)(scale * g);

        g = s = scale = 0.0;
        if (i < m && i != n - 1)
        {
            for (k = l; k < n; k++)
                scale += fabs((double)a[i][k]);
            if (scale)
            {
                for (k = l; k < n; k++)
                {
                    a[i][k] = (float)((double)a[i][k]/scale);
                    s += ((double)a[i][k] * (double)a[i][k]);
                }
                f = (double)a[i][l];
                g = -SIGN(sqrt(s), f);
                h = f * g - s;
                a[i][l] = (float)(f - g);
                for (k = l; k < n; k++)
                    rv1[k] = (double)a[i][k] / h;
                if (i != m - 1)
                {
                    for (j = l; j < m; j++)
                    {
                        for (s = 0.0, k = l; k < n; k++)
                            s += ((double)a[j][k] * (double)a[i][k]);
                        for (k = l; k < n; k++)
                            a[j][k] += (float)(s * rv1[k]);
                    }
                }
                for (k = l; k < n; k++)
                    a[i][k] = (float)((double)a[i][k]*scale);
            }
        }
        anorm = MAX(anorm, (fabs((double)w[i]) + fabs(rv1[i])));
    }

    for (i = n - 1; i >= 0; i--)
    {
        if (i < n - 1)
        {
            if (g)
            {
                for (j = l; j < n; j++)
                    v[j][i] = (float)(((double)a[i][j] / (double)a[i][l]) / g);

                for (j = l; j < n; j++)
                {
                    for (s = 0.0, k = l; k < n; k++)
                        s += ((double)a[i][k] * (double)v[k][j]);
                    for (k = l; k < n; k++)
                        v[k][j] += (float)(s * (double)v[k][i]);
                }
            }
            for (j = l; j < n; j++)
                v[i][j] = v[j][i] = 0.0;
        }
        v[i][i] = 1.0;
        g = rv1[i];
        l = i;
    }

    for (i = n - 1; i >= 0; i--)
    {
        l = i + 1;
        g = (double)w[i];
        if (i < n - 1)
            for (j = l; j < n; j++)
                a[i][j] = 0.0;
        if (g)
        {
            g = 1.0 / g;
            if (i != n - 1)
            {
                for (j = l; j < n; j++)
                {
                    for (s = 0.0, k = l; k < m; k++)
                        s += ((double)a[k][i] * (double)a[k][j]);
                    f = (s / (double)a[i][i]) * g;
                    for (k = i; k < m; k++)
                        a[k][j] += (float)(f * (double)a[k][i]);
                }
            }
            for (j = i; j < m; j++)
                a[j][i] = (float)((double)a[j][i]*g);
        }
        else
        {
            for (j = i; j < m; j++)
                a[j][i] = 0.0;
        }
        ++a[i][i];
    }

    for (k = n - 1; k >= 0; k--)
    {
        for (its = 0; its < 30; its++)
        {
            flag = 1;
            for (l = k; l >= 0; l--)
            {
                nm = l - 1;
                if (fabs(rv1[l]) + anorm == anorm)
                {
                    flag = 0;
                    break;
                }
                if (fabs((double)w[nm]) + anorm == anorm)
                    break;
            }
            if (flag) 
            {
                c = 0.0;
                s = 1.0;
                for (i = l; i <= k; i++) 
                {
                    f = s * rv1[i];
                    if (fabs(f) + anorm != anorm) 
                    {
                        g = (double)w[i];
                        h = PYTHAG(f, g);
                        w[i] = (float)h;
                        h = 1.0 / h;
                        c = g * h;
                        s = (- f * h);
                        for (j = 0; j < m; j++) 
                        {
                            y = (double)a[j][nm];
                            z = (double)a[j][i];
                            a[j][nm] = (float)(y * c + z * s);
                            a[j][i] = (float)(z * c - y * s);
                        }
                    }
                }
            }
            z = (double)w[k];
            if (l == k) 
            {
                if (z < 0.0)
                {
                    w[k] = (float)(-z);
                    for (j = 0; j < n; j++) 
                        v[j][k] = (-v[j][k]);
                }
                break;
            }
            if (its >= 30) {
                free((void*) rv1);
                fprintf(stderr, "No convergence after 30,000! iterations \n");
                return(0);
            }

            x = (double)w[l];
            nm = k - 1;
            y = (double)w[nm];
            g = rv1[nm];
            h = rv1[k];
            f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
            g = PYTHAG(f, 1.0);
            f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;

            c = s = 1.0;
            for (j = l; j <= nm; j++)
            {
                i = j + 1;
                g = rv1[i];
                y = (double)w[i];
                h = s * g;
                g = c * g;
                z = PYTHAG(f, h);
                rv1[j] = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g * c - x * s;
                h = y * s;
                y = y * c;
                for (jj = 0; jj < n; jj++) 
                {
                    x = (double)v[jj][j];
                    z = (double)v[jj][i];
                    v[jj][j] = (float)(x * c + z * s);
                    v[jj][i] = (float)(z * c - x * s);
                }
                z = PYTHAG(f, h);
                w[j] = (float)z;
                if (z) 
                {
                    z = 1.0 / z;
                    c = f * z;
                    s = h * z;
                }
                f = (c * g) + (s * y);
                x = (c * y) - (s * g);
                for (jj = 0; jj < m; jj++)
                {
                    y = (double)a[jj][j];
                    z = (double)a[jj][i];
                    a[jj][j] = (float)(y * c + z * s);
                    a[jj][i] = (float)(z * c - y * s);
                }
            }
            rv1[l] = 0.0;
            rv1[k] = f;
            w[k] = (float)x;
        }
    }
    free((void*) rv1);
    return(1);
}
//---------------------------------------------------------------------------
float *SolucaoSistemaSVD(float **A,int NumLinhasA,int NumColunasA,float *B)
{
   float **temp0,**temp2,**temp3;
   float *temp1;
   int i,j;

   temp0 = CriaMatriz(NumLinhasA,NumColunasA);
   for(i=0;i<NumLinhasA;i++)
      for(j=0;j<NumColunasA;j++)
         temp0[i][j] = A[i][j];

   temp1 = CriaVetor(NumColunasA);
   temp2 = CriaMatriz(NumColunasA,NumColunasA);
   DecomposicaoSVD(temp0,NumLinhasA,NumColunasA,temp1,temp2);
   for(i=0;i<NumColunasA;i++)
   {
      temp1[i]=(1/temp1[i]);
   }
   temp3 = MultiplicaVetorDiagonal(temp2,NumColunasA,NumColunasA,temp1,NumColunasA);

   temp1 = ApagaVetor(temp1);
   temp2 = ApagaMatriz(temp2,NumColunasA,NumColunasA);

   temp2 = CalculaTransposta(temp0,NumLinhasA,NumColunasA);
   temp0 = ApagaMatriz(temp0,NumLinhasA,NumColunasA);

   temp0 = MultiplicaMatriz(temp3,NumColunasA,NumColunasA,temp2,NumColunasA,NumLinhasA);

   temp2 = ApagaMatriz(temp2,NumColunasA,NumLinhasA);
   temp3 = ApagaMatriz(temp3,NumColunasA,NumColunasA);

   temp1 = MultiplicaVetor(temp0,NumColunasA,NumLinhasA,B);

   temp0 = ApagaMatriz(temp0,NumColunasA,NumLinhasA);

   return(temp1);
}
//---------------------------------------------------------------------------
float *SolucaoSistemaGAUSS(float **A,int NumLinhasA,int NumColunasA,float *B,int numElementosB)
{
   return(NULL);
}
//---------------------------------------------------------------------------
float **MultiplicaVetorDiagonal(float **A,int LinhasA,int ColunasA,float *B,int ElementosB)
{
   float **Resultado;

   Resultado = CriaMatriz(LinhasA,ElementosB);
   for(int i=0;i<LinhasA;i++)
      for(int j=0;j<ElementosB;j++)
      {
         Resultado[i][j]=0.0f;
         for(int k=0;k<ElementosB;k++)
            if(k==j)
               Resultado[i][j]+=A[i][k]*B[k];
      }
   return(Resultado);
}
//---------------------------------------------------------------------------
float *MultiplicaVetor(float **A,int LinhasA,int ColunasA,float *B)
{
   float *Resultado;

   int i,k;

   Resultado = CriaVetor(LinhasA);
   for(i=0;i<LinhasA;i++)
   {
      Resultado[i]=0.0f;
      for(k=0;k<ColunasA;k++)
         Resultado[i]+=A[i][k]*B[k];
   }
   return(Resultado);
}
//---------------------------------------------------------------------------
float  *CalculaAngulosEulerMatrizRotacao(float **Rw)
{
//Resultado[3] = { theta, *Tridente*, Phi}

   float *Resultado;
   float *Resultado2;

   Resultado  = CriaVetor(3);
   Resultado2 = CriaVetor(3);

   Resultado[0]  = -asin(Rw[2][0]);
   Resultado2[0] = M_PI + Resultado[0];
   if((Rw[2][0]!=1)&&(Rw[2][0]!=-1))
   {
      Resultado[1]  = atan2(Rw[2][1]/cos(Resultado[0]), Rw[2][2]/cos(Resultado[0]));
      Resultado[2]  = atan2(Rw[1][0]/cos(Resultado[0]), Rw[0][0]/cos(Resultado[0]));
      Resultado2[1] = atan2(Rw[2][1]/cos(Resultado2[0]),Rw[2][2]/cos(Resultado2[0]));
      Resultado2[2] = atan2(Rw[1][0]/cos(Resultado2[0]),Rw[0][0]/cos(Resultado2[0]));
   }
   Resultado2 = ApagaVetor(Resultado2);
   return(Resultado);
}
//---------------------------------------------------------------------------
float **CalculaInversaLU(float **Matriz,int NumLinhas)
{


   return NULL;
}
//---------------------------------------------------------------------------
float **CalculaInversaSVD(float **Matriz,int NumLinhas)
{
   float **Apoio,
         **Inversa,
         **A,
         **V,
         **Ut;
   float  *W;
   int     TamanhoA,
           i,
           j;

   TamanhoA = NumLinhas;

   A = CriaMatriz(TamanhoA,TamanhoA);
   V = CriaMatriz(TamanhoA,TamanhoA);
   W = CriaVetor (TamanhoA);

   for(i=0;i<TamanhoA;i++)
      for(j=0;j<TamanhoA;j++)
         A[i][j] = Matriz[i][j];

   DecomposicaoSVD(A,TamanhoA,TamanhoA,W,V);

   for(i=0;i<TamanhoA;i++)
      W[i] = (float)((double)1/W[i]);

   Apoio = MultiplicaVetorDiagonal(V,TamanhoA,TamanhoA,W,TamanhoA);

   V = ApagaMatriz(V,TamanhoA,TamanhoA);
   W = ApagaVetor(W);

   Ut = CalculaTransposta(A,TamanhoA,TamanhoA);
   A  = ApagaMatriz      (A,TamanhoA,TamanhoA);

   Inversa = MultiplicaMatriz(Apoio,TamanhoA,TamanhoA,Ut,TamanhoA,TamanhoA);

   Apoio = ApagaMatriz(Apoio,TamanhoA,TamanhoA);
   Ut    = ApagaMatriz(Ut   ,TamanhoA,TamanhoA);

   return(Inversa);
}
//---------------------------------------------------------------------------
